RECONSTRUCTION FROM LIMITED PROJECTIONS 


A Thesis Submitted 

in Partial Fulfilment of the Requirements 
for the Degree of 

MASTER OF TECHNOLOGY 


By 

AMITABHA SANYAL 


to the 

DEPARTMENT OF ELECTRICAL ENGINEERING 

INDIAN INSTITUTE OF TECHNOLOGY KANPUR 

AUGUST, 1981 



CcN ff. ■' 




08591 . 

I****** »*# «1f t«*»«>*4 trMVW 


ee-nsi-M- san-r&c: 



CERTIFICATE 


ii 



This is to certify that the thesis entitled 
"Reconstruction from Limited Projections' is a x'ecord 
of the work carried out under our supervision and 
that it has not been submitted elsewhere for a degree* 



(S.K. MULLICK) 


Professor 

Dept, of Electrical Engg.' 
Indian Institute of Tech:. 
KANPUR 


(r.k.s. rathore) 
Assistant Professor 
Department of Mathematics 
Indian Institute of Tech'. 
KANPUR 



ill 


ACKN OWLED GEMENT 


With a deep sense of gratitude , I wish to acknowledge 
my gratefulness to my thesis supervisors , Dr. S.K. Mullick 
and Dr. R.K.S. Rathore for their guidance, invaluable 
suggestions and encouragement throughout the course of this 

study . 


I also thank Shri Urinal Patra for his helps and 
suggestions. 

Finally, thanks are due to shri K.N. Tewari for his 
neat and efficient typing. 


si* 


\ J-v. 


AMITABHA SANYAD 1 


.1 



iv 


CONTENTS 


Page 


ABSTRACT 



Chapter 1 

THE PROBLEM 

1 

Chapter 2 

THE ALGORITHMS USED 

11 

Chapter 3 

SOFTWARE PROBLEMS 

32 

Chapter 4 

ENHANCEMENT 

37 

Chapter 5 

THE RESULTS 

44 

Chapter 6 

CONCLUSIONS 

56 

REFERENCES 


58 



CHAPTER 1 


THE_ PROBLEM 


The need for estimating the cross section of a object 
from its projections arises in a large number of areas. The 
most well known application is X-ray computerised tomography/ 
in which X-ray scans at different angles are combined to 
reconstruct a cross sectional image which represents the 
details of the scanned structure. The- performance of computer 
axial tomography in pinpointing tumors hemorrhages/ blood 
clots etc. was so dramatic as to revolutionize radiology [10*24-3 
A few specific applications are mentioned below. 

(1) Noninvasive angiograms; One of the traditional brain 
examinations / the angiogram, has the purpose of visualizing 
the arteries in the head, as well as certain tumors. The 
usual method of doing this is to inject iodine contrast dye 
into the cartoid artery in the neck. The dye goes directly 
to the head and then to the body'.' while the full amount is 
in the head an X-ray is taken. The direct arterial injection 
is painful and dangerous,, always calling for hospitalization/ 
and often producing grave complications. On the other hand, a 
slow intravenous injection of dye is considered relatively 
painless and harmless'.' But the concentration in this case 
is not sufficient for usual X-ray photography. However those 
concentrations are enough to be detectable in cross sectional 


reconstructions^ 
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(2) Even when a mass is known to be present in the breast# 
it remains a difficult problem to determine without, a biopsy 
whether the mass is cancer or a fibrocystic lesion. Since 
the fibrocystic lesion is more dense than the fatty breast 
tissue , 'While the cancer is still more dense# a reconstruction 
of the cross section of the breast can be achieved by simply 
examining the densitometer readings along a line through the 
lesion"." 

Apart from medical applications# there are other fields 
too,wher'' reconstruction tomography can be used extensively. 

In radio astronomy the two dimensional brightness distribution 
of a source can be estimated from fan beam telescopic scans. 

The problem can be extended to multidimensional fields too". 

In electron microscopy the three dimensional structure of 
biological specimens can be estimated from their two dimensional 
electron micrographs taken at different cilt angles. 

The cross sectional reconstruction can also be used 
to obtain three dimensional information by stacking successive 
transverse sections [30], Real ob rets exist in the three 
dimensional space R . However, a solution to the two dimensional 
problem provides immediately a solution to the three dimensional 
problem by recovering all two dimensional sections orthogonal to 
'a fixed line. The sections arc then stacked together to obtain 
a three dimensional reconstruction. At present almost all 
three dimensional reconstructions are being obtained in this 
way for reasons of computational simplicity. 
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Another potential application is in picture transmission 
[19] . in transmitting pictures # one can increase the rate of 
transmission of information faith only moderate loss of accuracy 
by transmitting only the projections instead of the actual 
picture. Consider a 256 x 256 picture. It has been seen that 
by reducing it to a coarser 32 x 32 picture and sending only 
32 x 16 projections# the picture car be adequately reconstructed. 
The compression of data in this case is 256 x 256:32 x 16 = 128 si 
Of course the accuracy can be increased by increasing the amount 
of information. But the transmitted data may be corrupted with 
additive noise and the reconstruction method should take into 
account this fact. 

In general# two types of geometrical arrangements are 
usually employed for producing the projections of a multi- 
dimensional field. These are based on illuminating the source 
by either parallel beams or fan beams of radiation (be they 
heat# sound# X*-ray or electromagnetic)';" Fig . 1 illustrates 
a parallel beam geometry. A set of collimated sources 
Illuminates the object and an arx~ay of sensors gather the 
projection data. By shifting the assembly of soux'ce and 
sensors through small angles, a large amount of projection 
data is gathered'.' In the case of X-ray tomography# the 
sources are X-ray guns and the sensors are low noise photo- 
detectors /which detect the amount of X-ray incident on them 
and the absorption in the path. 
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Parallel beam systems, specially in X-ray tomography, 
require very high precision moving parts to ensure correct 
registration of each projection at different angles. To 
minimise on rotating parts and to speed up the data collec- 
tion, fan beam systems have been recently introduced (Fig. 2) .These 
consist of a single rotating source which produces a fan beam 
of incident radiation. The radiation is collected after 
transmission through an object by a set of fixed sensors 
located on a circle. This is done for several angles and the 
projection data is gathered. In radio astronomy the scanning 
beam is essentially a fan beam. 

Mathematically , the reconstruction problem may be stated 
as follows [25] : 

(a) For parallel beam geometry (Fig.3)s Given a set of 
N projections of a two dimensional field or object [A(x,y)] 
at angles [©j,, k = 1,2,..., N] with respect to the: reference 
system of coordinates, i.e. 


B k (x k )dX k = )d X k dx k 

s 

where A represents the two dimensional field intensity at 
any point (x,y), 
y^ is orthogonal to x^, 

S represents the strip of integration of A(x,y) 
taken at the projection angle ©^ and location x^ and 
Xj. = xCosQj, + ySin© k , 
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* -xSinQ^ + yCos© k 


are the rotated Cartesian coordinates and Bj, ( x^ ) is a one- 
dimensional function with resoect to x, . The objective is 

K 

to obtain 'good' estimates of [A(x,y)] from observations of 
[ B k ( ^ ) ** 2Cj c / by appropriate signal processing* 


(b) Fan beam geometry (Fig.4)i The fan beam projection 
with the source at an angle 7 with respect to the reference 
system of coordin ates is given by 


where 


1 2 (0 b ) 

B r (e k ) = r k 

w 


A(x y y) dr^ 


and 


tanS = [.^°€aiJ[2±£2_tR ] 
* -xsinr + yCosr 

2 2 2 

r ]c = R + r - 2Rr Cos(r-0) 


R is the radius of the circle on which the source rotates 
and (r, 0) arc the polar- coordinates of any point (x,y).. 

0 k 2 * * * * 7 ^2^k^ are as s ^ own ti ' ie fi© ure ( Fig .4 ) . if 

the projections are taken at discrete angles [7=7^, k=l y 2, . . ,n] 
then [b, (0, )] represent the N fan beam projections of the 

■K /C 

object [A(x,y)]« 


One of the many points of interests is the extent to 
which the object A(x,y) can be recovered from its projections. 
The following theorems have been proved [24] in this connection. 



Fig, 4: ACQUISITION OF FAN BEAM PROJECTIONS, 
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Theorem 1, 1; An object is determined by any infinite set 
of projections'. 

JL\2 ; A finite set of projections tells nothing 
at all. 

Theorem 1.3 : Suppose given an infinitely differentiable 

object f and a finite number of directions. 

Then there is a new infinitely differentiable 
object f with exactly the same shape # exactly 
the same projections from these directions# and 
completely arbitrary on any compact set in the 
interior of the support of f . 

Theore m 1.4 t For almost any finite dimensional space F, the 
objects in F can be distinguished by a single 
projection from almost any direction'.' 

The second theorem is one of those mathematical ideals 
untainted by any possibility of practical application. Although 
a reconstruction method necessarily uses only a finite number 
of projections# it also necessarily produces only a finite 
dimensional space of possible reconstructions'. The question 
arises as to whether the projections serve to distinguish 
between the objects within this finite dimensional space. 
Theorem 1.4 says that the objects within this space can indeed 
be distinguished. 

So far we have formulated and dealt with the problem 
in its continuous form'. To enable the use of powerful digital 
processing techniques and the digital computer# the set of 
equations must be discretized. 
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For our purpose we have discretized the system so that 
for the ith projection, we get an equation of the form 


N 

E A . . x . *= B . 

j=l *3 J 1 


N is the total number of picture elements (pixels) in a 


NxN discretized picture. 

A . . = 1 or 0 depending on whether the pixel ’x,' is included 
3 3 

in the projection (ray) or not. 
is the sum of all pixels included in the particular 

projection. Thus we obtain a set of linear simultaneous 


eo nations 


Ax = B 

2 

where- A is a matrix of size MxN and £ is a vector of size 
Mxl; M being the total number of projections. 

The following points may be noted: 

(1) The matrix [a] is uite sparse, since from the geometry 
of projections, it is obvious that each projection 
includes vxry few of the pixels'.' 

(2) The size of the matrix [a] is enormous. In typical 

2 

applications N starts at 2500 and with little ambition 

6 9 

can easily reach 10 to 10 . The number of projection 

5 7 

elements ranges from 500 to 10 and 10 in corresponding 
cases. Thus the number of elements in the [a] matrix 
can range between 1.25 xlO to 10 • 
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(3) The ecru at ions are ordinarily highly underdetermined, 
i.e. there are less number of equations than there are 
unknowns".' 

(4) The rank of matrix [a] is unknown. 

(5) The matrix [a] is nonnegative since A.. > 1 or 0, 

(6) The projections are ordinarily nonnegative. 

(7) The pixel values x^ are assumed to be nonnegative 
so that one desires a solutiot for which > O'. 

(8) The errors in the data may cause the equations to be 
inconsistent',;" 

(9) Statistical effects of noise in the data should be 
analysed if possible. 

(10) The approximations which lead to the equations Ax = B 
introduce systematic errors which have to be analysed. 

In our case the approximations are: 

(a) All activity is assumed to be at the centre of 
the pixel. 

(b) A ray is assumed to be of zero width. No weight 
has therefore boon given either to the area of the 
pixel intersected by the ray or to the length of 
the ray that traverses the pixel. 

The projection data used for reconstruction is likely 
to get contaminated by noise. In the case of transmission of 
pictures, the data is likely to be corrupted by additive noise 
while it is being transmitted. In biomedical applications the 
X-ray fi 1m is affected by noise in the following manner [24]. 
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T'Jhen a photon from the X-ray beam passes through the 
object t T "0 things can happen. It can sail through completely 
undisturbed or it can interact 'with the atoms of the object. 

In the cas ' of an interaction the photon may be stopped which 
attenuates the X-ray beam. Unfortunately, however, the 
interactions are not this simple. Sometimes they involve 
not the stopping of X-ray photon but a change in its direction . 
Other times they involve the release of photons from the 
object atom. These new photons and the old ones with the 
wrong direction are called scatter.' There are many other 
sources of noise of which a f qvj ares finger marks of the 
technicians, streaks left by dirty rollers in the developing 
machine, incorrect developing temper a cure and uneven film 
emulsion". The noise may be large and comparable with the 
difference between the maximum and minimum pixel value. This 
may necessitate some preprocessing of the raw data. 

In all these cases the available data is a got of 
projections taken at different angles. Thus it is possible 
to Identify a gen or •-’1 problem, involving the reconstruction 
of an image from its projections, which is not restricted to 
the narrow? confines of any specific application. 

The aim of this thesis being the application and 
testing of a certain algorithm for reconstruction, we have 
not gone into an indepth study of the relative merits of the 
two typos of geometries for obtaining projections, but have 
restricted ourselves to the gathering of projection data by 
parallel beam geometry only. 
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T^e reconstruction problem has been solved by us by 
the use of a certain block iterative algorithm. The amount 
of data to be handled being enormous/ data structuring has 
been done for efficient computation and storage. The 
algorithm has been verified by taking a digitized picture 
and obtaining a linear set of simultaneous equations from 
it. The reconstructed picture has been visually compared 
with the original picture after printing both of them on a 
line printer, where the different grey levels have been 
simulated by overprinting. The mean square error has also 
been computed. 
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CHAPTER 2 


THE ALGORIT HMS U SED 


In Chapter 1 we discretized the problem to obtain 
a set of linear simultaneous equations of the form 

A x = B 

where A is a mxn matrix, x is a nxl vector and £ is a 
mxl vector'.' The methods for solving this system can be 
divided into the following categories s 

1) Direct, matrix inversion techniques? generalized inverse 
and pseudo inverse [ 8,12 , 27 ,2Sj . 

2 ) Summation, linear superposition or back projection 
techniques [30] . 

3) Algebraic reconstruction technique (ART) [5,10,11,15, 
16, 17, 30]. 

4) Simultaneous iterative reconstruction technique 
( c XRT ) [30,32]. 

5) Minimization of criterion functions [19]. 

6) Summation of filtered back projection, convolution 
technique [14, 17, 30] . 

7) Fourier reconstruction method [9,13,14,17,31]. 

1) Gen eral ized inverse and pseudo i nverse : 

A generalized inverse of a mxn real matrix A of 
rank k, is a convenient tool for handling numerically the 
system A x = B of m linear equations in n unknowns. Any 
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numerical procedure for solving such systems will have to 
distinguish four cases. Suppose that the rank of [A, B] is 

k' 

(1) k = k*. Then the equations have one or more solutions. 

i) If k = n < m, then there is a unique solution';" 
ii) If k < n, then the equations have an infinite 
number of solutions that can be written in the 

form x = x + u, i.e. u 6N(A) the null space of 

— - 1 

A, where N(A) is a vector space of dimension n-k. 

(2) k < k*. The equations are inconsistent* In this 
case it makes sense in many applications to look for a least 
square solution/ i.e. the solution that minimizes j| A x- B|| ^ * 
the size of the residual B - A x in the 2-norm. We again 
distinguish two cases? 

i) If k =n, there exists a unique least squares 
solution." 

ii) If k < n, the residual is minimized by a family 
of solutions = x Q + u, where x q is a particular least 
square solution and u € N(A) . 

The generalized inverse is motivated by the desire 
to cover all of those cases by one solution formula/ 

£ = x o + H' = 2 2' 5 being the generalized inverse. 

Def". 2,1 ; A g-inverse of an mxn matrix ^ is any matrix G 
such that 


A G A — A 
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T 


A symmetric g-lnverse satisfies A G A = A, (AG) = A G 


Def»~ 2V2 s A least squares solution of A x = B -is a solution 


that minimizes ||a x -B>| 
is one that minimizes 


A minimal least squares solution 


A X - Bt 2 and 


2 * 


T heorem 2,1 : A necessary and sufficient condition for x = G B 
to be a least squares solution of A x = 3 is that G be a 
symmetric g~ inverse.' 


Theorem 2.2 ; If G is a symmetric g~inverse, necessary and 
sufficient conditions for x = G B to be a minimal least squares 
solution of A x = B is that G also satisfy 

T 

G A G = G , (GA) - G A . 

Def« 2 .3 t The Moore-Penrose generalized inverse of A is a 
matrix G = A + that satisfies the four conditions 

A G A = A 
T 

(A G) = A G 
G A G = G 
( G A) — G A 

From the above definitions and theorems, it is clear 
that the Moore Penrose inverse may be a solution we are looking 
for'.' The methods of finding the Moore Penrose inverse are 
summarized below? 

i) If A is raxn (m ^ n) and it is known to be of full rank, 
the most straightforward method for computing A is via 
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A + = (A T A)" 1 


A (m > n) or 


T T • 
A (A A ) 


(m < n) 


This is a direct method'. The Greville algorithm is a 
recursive algoritlpm for computing the Moore Penrose inverse 
of a matrix. We present the algorithm without proof. 


(ii) Greville Algorithms Let A =(a^ ... a^ ... V 

where a^ denotes the kth column of the given matrix A and 

A* = % ••• V 


Consider A^ in the partitioned form 


a*. = V 


Compute 


and 


4 —k~l ^k 

% = ^k - 4-1 4 


If c-^ £ 0, then 



r _ 

-k -k' 


,-l 


t 

-k 


otherwise e^. = (1 + d^ d^_) ^ 


A* 

~ k-1 



is then given by 

r 4* 



} 

1 


^k 


»k -k 


To start the recursion/ take A^ = 0 if = 0 (null column 
vector); otherwise/ 






This algorithm requires only one decision in each recursion 
and is matrix inversion free'.' 
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2) Linear su perp osition or back projection 

wwwhumi f i mt 'i nmn ihl *»*usnm ******* tmMi&m+mm m*m m mimmm mm mmtmmmmrn 

The simplest and most rapid method of reconstructing 
a two dimensional distribution from multiple one dimensional 
projections is to merely project the views back to a common 
object region. Thus we describe back projection as 

x! = SB . 
i 3 

where the summation takes place over all those rays which 
have passed through x/. The total density or concentration 
for the section or array being reconstructed is estimated by 

T = SB . 

3 

where this summation is over all rays of a single projection 
angle. After back projecting, the total density T' for the 
array is 

N 2 2 

T* = S x! r N =n; the number of pixels. 

1 i 

A normalization factor is derived for reducing the 
value of each picture element so that the reconstructed array 
total density corresponds to the estimated total. Thus the 
corrected back projected image is 

x. = x f ( T/T 1 ) . 

1 l 

A more exact background correction involves modifying 
the values by subtracting from each pixel the mean density or 
concentration multiplied by the number of views minus one 

~ x [ ~ [T(no.of views - 1)]/N* 
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where N* is the number of sample points per projection. 


3 ) Algebraic reconstruction technique ; 

This simple algorithm consists of guessing a value of 
all the pixels and then compensating for the discrepancy 
between the measured ray sum Bj and calculated ray sum B j , 


x 


n +l _ n 


- x. 


B 


j 


If the measured ray sum and the calculated ray sum 
turn out to be the same then the guessed values are correct 
for the particular projection. However, for a different 
projection, there may be a large discrepancy. Thus the pixels 
ofthe old ray which also lie in the new ray willbe modified 
again. The above procedure is called multiplicative ART. 


There is another form of ART in which the discrepancy 
between the measured and calculated ray sum is corrected by 
adding to each pixel the difference between the two divided 
by the number of pixels in the ray. 


x 


n +l r n B j “ B i pi 

= max x. + *• , CJ 

i *- "i 


j 


Nj is the number of pixels in the ray which passes through 
the pixel xd. This is called additive ART. 


There are a few modifications of ART. One consists 
of setting to zero those values that were zero, since they 
corresponded to a ray sum that was zero. This is a boundary 
condition for this iterative technique. Another form of ART 
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takes into account noise and has been used effectively in 
transmission studies and simulations with noise. 


It may be noted that the following unconstrained 
version of ART is the same as the iterative matrix inversion 
method by Kacmarz. 


x P+1 = x P 


N" 


Cj + [A. .(B.- E A. . x . )]/N. 

J 0 AJ i i=1 ij J J i 
N 2 2 1 

~ * A ±1 - / 1 = 1/2/ .. . m. 
j=l J 

4) Simultaneous iterative reconstruction technique; 


where 


N. 

l 


The simultaneous iterative reconstruction technique 
differs from ART in the fact that at each iteration the pixels 
are modified by using data from all projections simultaneously. 
Thus 


„ , ■, „ SB. SB • 

x! = max[x. +. - — i-/0] 

1 1 SL j SN. 

J J 

Here the summations ar 5 done over those rays which intersect 
the pixel".' Lj is the length of the ray jV B^ is the measured 
project density of ray j, . is the calculated projected 
density of ray j,and Nj is the number of pixels in ray j. 

A multiplicative algorithm has also been given , which 


is 


,n+l 

x i 


max[ 


SB. SN. 

I 

Sl. SB'. 
J J 


n 

v °] 


After each iteration the total array is normalised such that 
for all i x? +1 = (T/T 1 ) x| n+ ^/ where T and T* have been 


defined before". 
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5 ) Minimization of criterion functions s 

An example of this method is the algorithm suggested 
by Kashyap and Mittal in which sum of two quadratic functions 
and has been minimized subject to the condition A x = B. 
is the nonuniformity function 


1 _ , - *2 it 

2 f < x i - x i> = 2 .2 £ 2 - 

i consists of all cells not on the boundary/ x^ is the 
average of eight neighbours of cell i and c is a positive 

definite matrix. J 0 is the loss function 

2 Z 

. N , o 

■J S (x ± - a) z 
i=l 

2 

where a is the mean of intensities of the N cells. 

, T IT 1 T ' 1 . t 2 2 

J 1 + J 2 = 2 - c x + ^ ^ x ~ 2 N a 


_ - 1 at 2 2 
= J 3 - 2 N 3 


If we assume that a is known, then the problem is to minimize 
j o subject to A 25 = B» For this a function J(x) is chosen 


J(x) 



T IT 

x C x + 2 ~ 



2 


in which Y can be chosen experimentally and p is a large 
positive integer’.' J(.x) is minimized by a recursive method/ 
the steepest descent method. 


x(t-fl) = x(t) - (pfij/djc) 


x = x(t ) 


= x(t) - p[x(t) + Y c x(t)- P (y-B 5 j(t) )] 

t = 0/1/2... 


p is a scalar gain which can be chosen to ensure proper 
rate of convergence';" 
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^ £2H£ig j£. Reconstruction Method 

To describe the Fourier reconstruction method/ we 
shall introduce a new set of notations. Let f(X/y) or 
f(t/ 0 ) be the intensity distribution of the object and 
g(lz ©) be the projection at position 1 at angle 9. 

Then f(x,y) or f( r/ 0 ) ca n be expressed in terms of g(l/9) 
by the following Fourier transform relationships. 

co 

F(R# 6 ) e. / g( 1/0) exp(~2TriRl)dl (2.1) 

-QD 

2n cd 

f(r # 0) ~ j j f(R/ 9) exp[ 27iiRr c os ( £5-9 ) ] RdRdG 

° ° ( 2 . 2 ) 

The Fourier reconstruction method consists of either 
directly evaluating these integrals or using one of the 
methods described below s 

The basic equations can be rewritten in a form not 
involving Fourier transforms/ but containing only integrals 
of functions defined in the real space of observation';' We 
can write eqn.( 2 „ 2 ) as 

71 OO 

f (r/0) = / f j R j F(R,9)exp[27TiRr Cos(0-9)]dRd@ 


o — QO , 

(2.3) 

If we define. 



> , oo 



g* (1/9) ~ . J* j 

r| F(R/0) exp(2TtiRl)dR 

(2.4) 

-oo 



then from equation ( 2 '.' 3 ), 

71 

f(r # 0) ss f g • [ rCos ( 0-9 ) / ©] d0 
o 


(2.5) 
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Fourier inverting equation (2.1) we get 

0D 

g(l,e) = I 

-CD 

Comparin g equations (2.4) and (2.5) we obtain 

F.T. of g' (1/©) = F.T . of g(l,©) x F.T. of q(l) (2.7) 

where i is the Fourier transform of q(l) or 
co 

J R } = / q(l) exp(27TiRl)dl 

-CD 

It may be noted that back projecting g*(l/©) will give us 
the true image ( eqn .2.6)'. There are two methods' to evaluate 

g*(l,e) * 

a) Convolution Method: 

From the convolution theorem for the inverse of the 
product of Fourier transforms and equation (2.7) ,we can write 

oo 

g' (1/0) = / g<l ,©) qd-l-^ dl 1 (2.8) 

-co 1 

For calculating g'(l,©) we require q(l) explicitly. 

Now, 

oo 

q(l) = / t R | exp(-2 7r i R 1) dR 

-OO 

This integral cannot be evaluated between the limits -~oo 
to cq as the integrand diverges* - So to overcome this 

i 

difficulty, we replace -oo and -foo by -A/2 and +A/2 where 

A is a finite number. 

A/2 

q (1) = / (r| exp(2xr i R l)dR 

Jr\ & 1 

-A/2 


F(R, ©) exp[27TiRl] dR 


( 2 ,. 6 ) 
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2 

Thus we get, q^(na) = 1/4 a for . n = 0 

= ~l/n 2 n 2 a 2 for n odd 

= 0 for n even, 

q(na) being a set of points at equal spacing a. 

Equation (2,8) can be expressed in the form of an infinite 
sum 

g'(na,9) = a E g(ma, 9) q[(m-n)a] 
m=-oo 

= g(na, 0)/4a 2 - (1 /tt 2 ) E g[ (n+p)a,9]/p 2 

p odd 

The true image can now be generated by back projecting 
g*(na, 9). 

Filtered back projections 

The second method, called filtered back projection, 
starts with the following equation, 

F.T. of g* (1,9) = F,T. of g(l, 9) x !R| 

Sometimes instead of j R [ we have w(r). |Rj, where w(r) is 
a v/indow function"! The R.H.S.is now inverse Fourier trans- 
formed to give g'(l,9). This is now back projected to yield 
the true image! 

It may be noted that both these methods are essentially 
the same. The first method carries out the operations in the 
real space of observation whereas the second method operates 
in the frequency space"! 



The Method Used 
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This algorithm proposed by Rathore , is a general 
stationary least square method to solve a rectangular set 
of equations 7 This is an iterative method and its 
convergence under very general conditions has been estab- 
lished.' Even for inconsistent systems, the method provides 
a least square solution. One good feature of this method 
is its guaranteed convergence. The method is now described. 


Consider the ecu at ion A x = B, An asumed vector 

o o 

x will give rise to a residual error vector £ given by 


o _ 1 , o 
r_ = B - A x 

f «*0 — 


Thus r° = A[2£ - £°] 


= h. 


where e is the error vector'. 
- o 


If we can somehow calculate the error vector approxi- 
mately, we get a improved solution 


2£ 


o o 
2 £ + 


o o 

where e, is an approximate solution of r 




The iterative method' may be continued and the process 
has been shown to converge. The method of residual projections 
comp.utes the error vector by projecting the residual at every 
stage on the columns of matrix A heuristic reason for 
doing so is as follows: 
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The solution ^ of the problem satisfies 


B = A X 


n 

B = £ 

j=l 



x . 


3 


where a^ are the column vectors of A. The values x^ may 
be assumed to be the contributions of each column vector 
towards the making of B'.‘ The solution 2 is a set of contri- 
butions for which the residual £-A x is minimized in a 
least square sense. When the residual is orthogonalized 
with respect to the first column vector of A, it means that 
the component of the first column vector along the residual 
is zero. This is done for every column of A. The vector x 
for which the residual has no component along any of the 
columns of A is the required solution. 


The residual projection algorithm consists of the 
following steps: 

a) Start with any arbitrary initial vector and compute the 
corresponding residual. 

b) Project the residual computed in (a) on the first column 
of matrix A and compute its component along first column. 

c) Alter the first entry in the initial vector by adding 

the scalar component of the residual computed in (b) to it. 

d) Remove the vector component of the residual along the 
first column to obtain new residual. 

e) Project the residual obtained in (d) on the second column 
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of A ana update the second entry in the solution vector 
by adding the scalar component of the residual along the 
second column to it. Subtract the vector component of 
the residual along the second column : .from the residual 
obtained in (b ) 

f) Continue the process till the last column is finished. 

g) Go to (b) and repeat the process through (f) till 
satisfied. 

These steps can be mathematically expressed as 
follows : 

i) Start with an initial vector x° 

... o *• o 

11 ) Compute the residual r Q = B - A x 

Then for each k = 0,1/2,... till satisfied do, 
iii) For each j = 1,2,... N do 

£j = £jl x a J/Uj, aj] 

and x j +1 = x j + [ — j ^ ^ j * — 

k keeps track of the number of iterations and j keeps track 
of the column on which the residual is being projected. 

The following theorem has been proved. 

Theorem 2.4 : The vector sequences [ry j i = O ...00 and 

[x (i) ] i = 0 ... 00 , generated by the method of residual 

projections converge for arbitrary A/B and x(°) with 

x (ao) = ij.m r^ = P, B and x^ 00 ^ = lim x^^—Px^°^+G B. 
i-> on ~ ~ i— * 00 “ 
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Moreover, x^°°^ minimizes B-A x and G is a reflesive least 
squares g inverse of A. £ is the projection onto kernel 
of A along image of G A and is the orthog ;nal projection 
onto kernel of A. 

Residual projection with averaging 

Since the system of equations taken*' for reconstruc- 
tion was almost always underdetermined, the least square 
solution provided by the residual projection algorithm did 
not always coincide with the desired solution";’ The restoration 
often resulted in ‘ghosts'; pictures entirely different from 
the actual picture. To resolve this problem, additional 
information had to be provided in the form of physical 
constraints'." 

Inspired by the idea of Kashyap and Mittal [19] 
the idea of continuity was applied to our method".' Any 
real picture is continuous and the value of each pixel is 
nearly equal to the average of the pixels surrounding it. But 
this property, as applied by Kashyap and Mittal, gave rise to 
large matrices whose elements were other than one and zero. 

Thus efficient computation and storage was impossible. Each 
iteration took 30 to 40 seconds by their method. Their idea 
was modified to suit our purpose";" 

i 

Instead of minimizing the quantity 

1 { * - ,2 
j = -i S '•x. - x. ) 

z 16s 1 
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where s is the set of pixels not on the boundary and 
is the ' average of the neighbours around x ^ , we tried to 
force each pixel (not on the boundary) to a value which was 
the average of all the eight neighbours around it. This 
was done after each iteration and the values of x thus 
computed were used to find the new residual B~a It may 

be noted that the implementation of this scheme is extremely 
simple, needs almost no storage and is extremely fast.' 

Thus after each iteration we have 

x. = ~ E Neighbours of x . , 

i o x 

where x^ is not a boundary pixel. Thus we get the following 
matrix equation 

2 = C x 

2 2 2 

where £ is a N x N matrix and N is the total number of 
pixels. The elements of C are as follows. 

If is a boundary cell then 

C . . = 0 for i j 

ij 

= 1 1 for i = j 

If x.£ is not a boundary cell then 

C . = 1/8 for j=i— N-l* i-N/i-N+l/i-l/i+l/i+N-l* 
ij 

i-tN, i+N+1 

- o for all other values of j.‘ 

This operation has been used in conjunction with the 
residual projection method. The overall method converges 
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if the .iteration matrix of the compound operation has a 
spectral radius less than 1.' 

Residual Projection with Median Filtering 

A variation of the strategy mentioned earlier is to ' 
do median filtering at each iteration. In recent years 
median filters have come into use for noise reduction pur- 
.pose® in signal and picture processing. The main merit of 
a moving median is that it preserves edges, in constrast to 
a moving average which smooths edges. 

The median of n numbers Xj, , . » x^ is, for n odd, the 

middle number in size and we denote it by Median (x, . . .x ) , 

In 

For example the median of 0, 3, 4, 2, 7 is 3. 

A two dimensional median filter with filter window 

i 2 

A on a picture[x ±j ./ (i,j) 6 z ] i s defined through 

y ij = Me J ian x ij = Median (x ±+r ^ +s ; (r,s)€A), (i,j)e z 2 

Various forms of filter windows A can be used, e.g. line 
segments, squares, circles and crosses. For such symmetric 
filter windows, the median filter is expected to preserve 
edges'," The median filter also removes impulse and salt and 
pepper noise if the rate of error in the noisy picture is low. 

In our experiments, a 3x3 rectangular window was used. 
Since median filtering preserves edges, so should iterations 
of median filtering. Median filtering was done at each 
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iteration of the residual projection method to arrive at 
the' exhibited result. 

Computational Efforts 

The computational efforts required to carry out 
the different methods of reconstruction have not yet been 
calculated in terms of the picture size, the number of 

v 

equations and similar parameters'. An effort has been made 

s 

by us in this direction. 

Let the picture size be NxN. Let m be the number of 
angles and k be the number of rays per angle, k is usually 
taken to be equal to NV The total number of equations 

is mk. Since the matrix A is sparse, it is essential that 
the computations be carried out with the nonzero entries 
only'.' A software implemei tation of such a strategy is 
presented in the next chapter 7 

An important problem is to first relate the number 
of nonzero entries *p* of the A matrix to the parameters 
N, k and nu When the rays are parallel to the sides of the 
picture we expect that each pixel will be intersected by 

one ray in a particular angle. For such angles the number 

1 2 
of non-zero entries per angle is N . 

For other angles some pixels will be left out. If 
the ratio of the area of the picture not intersected by any 
ray to the total area of the picture is a, then the number 
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• . 2 - 
of pixels not intersected by any ray will roughly be <XN ! 

Simple trigonometric calculations will show that this 

ratio <2^ for any angle . $ is 


a. 


tan p 


[l— tan J3 + — ~ ] 2 

Cosf? s inj3 


The value of a is maximum when the angle • (5 is 4 5°', 

and this value is 0.086! Therefore p will be of the order 
2 - 

N m'.' This has been verified in our experiments. With a 

3 

32 x 32 picture and 17 angles we had 16x10 nonzero entries'.' 
Similarly with a 16x16 picture and 9 angles the number of 

3 

nonzero elements were 3.5x10 . 


Having given a rough estimate of p we show the 
computational effort of each method in terms of p'. 

1) ART : This method requires P multiplications and p 

additions for calculating the ray sums and p 
multiplications for modifying the pixel values'! 

This is for a single iteration"! 

2) SIRT : The computational effort for SIRT is difficult 

to calculate! We give the figures for the worst 

case when all the rows of the A matrix contain 

2 

either all zeroes or ones! Then pN multiplications 
2 

and pN additions are required to calculate the 

2 

sum of ray sums and N multiplications are needed 
to modify the pixel values! This is also for 
a single iteration. 
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3) Back Projection ; This method requires p multiplications 

and p additions for projection. In addition p 
multiplications and p additions more are required to 
calculate the normalization f actor V Thus 4p 
operations are required in total." 

4) Convolution Method % Here the computational requirement 

2 2 

is k m multiplications and k m additions for convo- 
lution and 4p operations for back projection". 

5) Filtered Back Projection ; The filtered back projection 

Ion 

method requires — log^ ^ complex multiplications and 

km log^ km complex additions for computing the FFT of 

]gn 

projections, km multiplications for filtering, log^ km 
complex multiplications and km log 2 Jen complex additions 
for computing I FFT and 4p operations for back projection. 

6) Residual Projection Method 2 This method requires p 

multiplications and p additions for projecting the residual 

on the columns of A matrix, p additions for modifying the 

2 

residuals and N multiplications for modifying the pixel 
values. These are the values for a single iteration".' 

Reconstruction with noise ; Since the process of averaging 
reduces noise, the residual projection method along with 
averaging should be ableto reconstruct pictures in the 
presence of noise. 
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We expreimcnted by adding uniform random noise 
between 0 and 1 to a picture whose maximum intensity was 5. 
Projections were obtained from the resulting picture and 
an attempt was made to reconstruct the original picture 
from these projections'. 
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CHAPTER 3 


SOFTW i RE PROBLEM S 


The major hurdles in a reconstruction problem are 
computation and storage. Even considering a 64 x 64 picture, 
the number of unknowns is 4096 and the weighting matrix A is 
4096 x 4096, if the system is to bb of exact deterrninacy. 
Direct inversion techniaues are therefore out of question. 

The method' must use the total available data in parts. Also, 
the sparsity of the A matrix has to be made use of. 

As an example consider the case of a 16 x 16 picture 
solved by the method of residual projections. The number of 
angles wore 9 and the number of rays per angle were 16, which 
gave a total of 144 equations. Therefore a A matrix whose size 
256 x 144 i.e. 36864 elements had to be stored. It can easily 
be seen that increasing the picture size further can pose 
enormous problems'.' 

A way of resolving this problem is to either produce 
or read each column separately into the core, use it for 
computation and then call the next column. The residual 
projection method is particularly suited for such a scheme. 

But then this means that there are 256 read statement 
executions or subroutine calls for each iteration. Since a 
disk read or a subroutine call takes much more time than a 
basic operation like multiplication, the whole process takes 
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enormous time. Thus we have saved in memory space at the 
expense of time. 

If we can load only the nonzero entries of the A matrix 
into the core and devise a method in which the only operations 
are performed with these entries/ then a saving in time and 
memory space can be achieved. Hew this is exactly done is 
shown below. 

v Advantage is taken of the binary nature of the A matrix, 
The elements of the A matrix are numbered thus t 

F 9 ! 

j 1 M+l . . . (N ~1)m+1 

| 2 M+2 

3 

4 

• # «* 

# * • 

j 

M 2M N 2 M I 

Now the A latrix is produced columnwise and if/ say, 
the 5th element is the first nonzero element (1) then the 
first entry of a one dimensional matrix IA i.e. IA(1) = 5. 

If the second nonzero entry is the M+lth element then 
IA(2) = M+l. Therefore if a A matrix is of the form 

10 0 o’ 

0 10 0 

0 110 

0 0 0 0 

1 0 0 0 ^ 
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T h en -it will be stored in the IA vector in the form 

IA (1,2,. 1/4/5) = 1,5,7,8,13. Storing in this manner., 

we have round that the total size of the IA vector in the 

case of a 16 x 16 x-icture has been approximately equal to 

3 550V A reduction of 1/10 has been achieved. A slight 

modification has to be done to the IA vector to find out 

whether the last entry of the A matrix has been reached. One 

2 

more entry whose value is more than N M, is added to the 
vector IA. In the previous example, the IA vector will be 
IA(1,2,3,4 ,5,6) = 1/5,7,8,13,21. 


To show how such a storage mechanism could be used 
in our algorithm, we have shown the steps for one iteration 
in details. 


The steps, expressed mathematically ,^re 


For j = 1 to N do 


1 ) 6 . = nT r , 

j ~j ^j-1 


whore Dj is a column vector of matrix D 

h isascalar - 

2) xj - x. - 

3) £ 1 -i-i 


and' is oqu^l to Pj^j 


6 . a, 
J -J 


In algorithmic language and using our storage strategy we 
write the same thing as follows : 
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1) m' = 1 

2) Do for j = 1 to the steps 2 through 17. 

3 ) n ' = m' 

4) 6j = ° 

5) j* = j x M 

6> Compare IA(m') with j* 

If IA(m') <: j' go to step 7 
If TA(m') > j' go to step 11 


7) 

i ’ = IA(m' ) - ( j — 1 ) M 


8) 

6 j = 6 j + r i‘ p j 


9 ) 

m" = m* + 1 


10) 

Go to step 6 


11) 

x j = h + h 


12) 

Compare XA(n') with j* 



if IA(n*) < j 1 go to step 

13 


if IA(n‘) > j’ go to step 

17 


13) 

i* 

= !A(n* ) - ( 

j“l) M 

14) 

r i- 

= r. , - 6 . 
i' J 


16) 

n * 

= n' + 1 


17) 

Continue 



Similar strategies have been aborted for carrying out 
all the operations prior to starting the iterations. It 
should be noted that all matrix operations like matrix 
multiplication etc. must essentially be done columnwise. 
This strategy reduces the number of computations from 
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2 2 

0(N M) to 0(kN M) where k is a fraction. (In our example, 
k = 0.1 for the 16 xl6 picture and 0.01 for the 64 x 64 one) . 
The value of k depends on the picture and the angles of 
projection. Further reduction may be made if a symmetry 
is noted in the position of the entries of the A matrix. For 
the pictures and the particular projection angles that we 
havo taken no such symmetry was observed. 

It may be noted that even if attenuation correction 
schemes are used, resulting in the entries of the weighting 
matrix [a] being other than 1 or 0, the order of computations 
or storage will not change. An additional matrix VALIA whose 
size is the same as IA will have to be used to store the values 
of the different nonzero entries of the A matrix. The matrix 
IA ill still indicate the position of nonzero entries only. 
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CHAPTER 4 
ENHANCEMENT 


ilien we solve an undsrdetermined system using the 
iterative .algorithms, the reconstruction in some cases may 
turn out to bo poor. Furthermore/ when a physical constraint 
like continuity is applied, the picture may be blurred and the 
edges may loose their sharpness. This is all the more true 
in our case since each pixel is the average of the points 
around it. So for visual satisfaction some operations must 
be done on the reconstructed picture. 

Several methods have been tried by us to varying 
degrees of success. Visual examination revealed the following 
errors in our reconstructed picture; 

1) The edges were not sharp due to the process of averaging. 

2) Th rest of the picture was riddled with salt and pepper 
noise (isolated dark points in light regions and vice- 

versa) . 

Thus the enhancement method had to do the following: 

1) Smooth the picture without blurring the edges. 

2) Sharpen '.ho edges'.' 

Removal o f salt and pepper nolse f3ls?o smooth a picture 
without blurring it, one can use nonlinear operations that 
involve averaging. The simplest class of such operations 
combines averaging with thresholding . One can clean up suit 
and pepper noise by making a point 'whiter* if the average 
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gray level in some neighbourhood of it exceeds its gray 
level by more than some threshold. The procedure is similar 

for mebing a point blacker'.' 


T .‘ r e have tried a simple implementation of this scheme. 
For a particular pixel x. / let the average of its eight 
neighbours bo x . . Nov/, if the quantity x x . is higher 

j j j 

than 1/4 th of the maximum possible value, then x^ is 
replaced by Xj? otherwise Xj is left unchanged. This process 
does not blur the picture but tends to smooth it out. 


Sharynlng the <3dgg[ 3];Since integration (or averaging) blurs 
a picture, a natural approach to deblurring is to perform 
some sort of differentiation operation. If deblurring is 
ronui rod only in some particular direction, one can simply 
take a directional derivative , For most purposes, however, 
one wants to dob] ur in every direction; this can foe accom- 
plished by taking the derivative in the gradient direction, 
i.c. the direction in ■ hich the grey levels change fastest 
at each point. For a well behaved function, this maximal 
directional derivative is equal to the square root of the 
sum of squares of the derivatives in any pair of orthogonal 
directions, o.g. [(Sf/6x) + (6f/6y) ] • For a digital 

picture one would approximate the derivatives by differences. 


e.g, 


[ (a 


i+l,j 


- a ) ' 


+ (a i,j+3- 


2 1 /2 

Vi 
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Similar results are obtained by using a sum of absolute 
values of. differences rather than the square root of the sum 
of their squares. Typical digital gradients are 


a> + Vl,] + - < a --- - + a -- .. . , ) 


► j“^l 1+1 /j~l 1+1/ j 1+1 /j+l" 


+ f ( 0 l - . - + 3l . , +a ) 

2-1/ j-1 1/ j-1 i+l/j— 1 


+ a i/j+l + a i+l,j+l^ } 


1/2 


b) (a 


i-l/j-l + a i-l,j + a i-l / j+l^ ^ a i+l,j-l + a i+l, j +a i+l, j+1^ 


"t (a, , + a. . - + a. . .) 

i-l/j-1 i,j-l l+ljj-l 


(a. + a. , , + a. . . _) 

i~l r j+l i/j+1 l+l t j +1 


Another useful combination of derivatives is the 


? q 2 2 

La pi act an (d"f/6x ) + (6 f/Sy ). For a digital picture the 


aru'tloqr’ts oppression is 


C (a t,j " Vl.j’ " (a i+l,j ' a ij^ 


■ij "i.J-l' "X J+l ~ h j ^ 


+ On - a i .i-i 5 1 (a i 


« 4a 


. . - (a. , , + a + a . , + a. , ) 

Ij i~l i J 2+1 / J 2/J* ,,, l 2 / J +1 


Using the Laplacian a filter was made, which was 
supnosed to enhance the edges. The f i 1 ter performed the 

following operation: 


6 2 . 6 2 


y(x 1 , x 2 ) « x(x 1 , x 2 ) + 0.2(~-j~ + ~~^~2 )x(x x > x 2 ). 


1 / 2 


6x- 6 x 2 
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•^P’py.PP - . P-- Since we worked with an 
und*.jid iiTuned system, the solution was not exact and some of 
the piA^ls which nad lower values in the original picture, 
oven turned negative. These pixels usually formed the back- 


ground and to clean it up we replaced all negative entries 
by zeroes . A variation of this is to replace negative entries 

by zeroes at each iteration. 


Kdg o .Pppy-py.fo ? smoothing by using masks [2] There have been 
many no ore on the subject of smoothing a digital image. A 
basic trouble with smoothing is that it tends to blur any 
sharp edges which happen to be present. This edge preserving 
smoothing method attempts to resolve the conflict between 
n o : se elimination and edge degradation. It looks for the 
most h' O’, m- . -no.-iur. region around each point in a picture and 
then gives each point the average gray level of the selected 
noiqhbc. . '/hood area. Noise in a picture is removed by repeated 
use of thin nv thoa, while the edges remain sharp. 

The procedure of the edge preserving smoothing is as 

follows : 

1) Rotate an elongated bar mask around a point (x,y). 

2) Pot act the position of the mask whose variance of the 

pray 3 cvel is minimum'. 

3) Give the average gray level of the mask at the selected 

posit. ’“on to the point (x,y). 

4) Apply steps 1-3 to all points in the picture. 
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5) Iterate the above method till the gray levels of 
the picture do not change. 

dtr to remove the noise without blurring a sharp 
edge, averaging must not be applied to the area which contains 
the o "■.go . Ii an area contains a sharp edge, the variance of 
the gray .1 cv--\l in that area becomes large. Therefore we can 
use the variance as the measure of nonhomogeneity of an area. 


c uppose that a picture has two regions R^ and R 2 , whose 

2 2 

mean and variances are (0, c^) and (m, o 2 ) respectively. Let 

a point (x,y) belong to R^. If (x,y) is located in the 

central portion of , its gray level approaches the gray 

level of after sev-ral iteration. Hov.ever, if (x,y) is 

near th*"* boundary of there exists two kinds of neighbourhoods, 

one which completely includes and the other which contains 

2 

part of R ^ and R 2 - The variance of the former is about . 

The variance of the latter can be calculated as follows. Let 
J ‘l ^*2 c '‘ :no ‘ 0 the number of points in R^ and R 2 respectively, 

which are contained in this neighbourhood. Then the variance 


N, 


M 


*** 

[S C 

1=1 


1 ^2 \ 2 , N 
x i ~ r m) + j=1 


2 2 
S^(x. 


N„ 

N 


m) Z ] 


1 2 

whore N ” N 1 + N 2 and x ± and x. denote the gray levels 
of the points belonging to R 1 and R 2 respectively. Expanding 

the R.H*s.,we have 


1 

N 


> 


2 

l 0 ! 


2 

N 2 0 2 




2-i 

m J 


N 
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N 1 1 2 2 

where we used S (x.) = N cr , 

1=1 1 1 1 

N 

- 1 

Ex, = 0 , 

i=l 

2 2 

3 Z =1 (x j * m) = N 2 °2 ' 

N 

2 

and E ($ . ) = N 0 m. 

J=1 J 2 

If a 2 < a 2 , that is a 2 + (Nj/N)m 2 < o 2 the 
neighbourhood containing both parts of R^ and are selected. 
Then the boundary between R^ and R^ is blurred by this opera- 
tion. But in most pictures, it is reasonable to assume 

2 2 2 2 

a l = °2 and 0 > a l* Then the correct neighbourhood is 

selected. 

For actual realization, we have used the nine masks 
shown in the figure (Fig. 3). The variances of these nine 
masks are compared with each other and the average gray level 
of the least variance mask is given to the point (x,y). 

The fact that the above method not only smoothens but 
also sharpens a blur-red edge can be understood by considering 
the following one dimensional example. The following is a 
one dimensional digital blurred edge 
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Gray level 0 

Mean 0 

Variance - 

(x 3) 0 


o. 

0 

0 

1 

3 

5 

8 

8 

8 

0 

0 

1/3 

4/3 

3 

16/3 

7 

8 

8 

0 

0 

4/9 

14/3 

8 

38/3 

6 

0 

0 


The mean and variance have been calculated using the 
1x3 mask centered at each point. If we select the minimum 
variance mask and give the average gray level of the selected 
mask to each point we get the followings 

0.000017888 

where the average gray levels 'have been rounded off to the 
nearest integer. This ^repeated once again,gives a completely 
sharpened edge.' 

00000088 8 8 

The fluctuation of the gray levels is gradually reduced 
by several iterations of smoothing. Once a point has a neigh- 
bourhood of a constant gray level, its gray level is never 
changed, by smoothing. Therefore, the number of points whose 
gray levels are changed by smoothing will gradually decrease 
to zero. The number of iterations needed for convergence 
depends on the amplitude of noise fluctuation and the shape of 
regions in the picture'.' 
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CHAPTER 5 

THE RESULTS 


The results of the various simulations carried out 
are presented in this chapter. The programs were run on the 
DEC 1090 computer and the pictures printed on the line 
printer. The gray levels were simulated by overprinting 
different sets of characters. A mean square error 
criterion was also computed for each of the reconstructed 
pictures. If x^j and y^ . are the original and reconstructed 
pictures of size NxN.then the mean square error MNSQR is 
given by 

, N N ~ 

MNSQR = -w'/'f S S (X. . - y .,) ] 

N 1=1 j=l 1J 13 

The error criterion MNSQR gives a reasonable estimate 
of the accuracy of reconstruction'." 

The residual projection algorithm guarantees a decrease 
of the residual I f B-Ax | j 2 » This has been verified and the 
residual has been tabulated in Table 1. This is for the 
32 x32 picture of Fig.l. 

Various initial conditions were tried".' The one that 
we found most satisfactory was the following. For each pixel 
the number of ones in the corresponding column of the A matrix 
is counted! T his number multiplied by 0.1 is the starting 
point for the particular pixel. 
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TABLE 1 


Iteration No '.* 

Residual 

Iteration No. 

Residual 

% 

1 

546.01 

16 

15.85 

2 

3 52.14 

17 

12.56 

3 

256.14 

18 

10.94 

4 

196 *.4 7 

19 

8.92 

5 

166.92 

20 

7.41 

6 

128.82 

21 

6.68 

7 

107.83 

22 

5.98 

8 

85.27 

23 

5.57 

9 

69.16 

24 

5.20 

10 

59.10 

25 

4.83 

11 

45.72 

26 

4.60 

12 

38.15 

27 

4.41 

13 

31.16 

28 

4.27 

14 

23 .33 

29 

4.17 

15 

19.65 

30 

4710 


Figs. 1 to 5 are the original pictures used for 
reconstruction';* Fig*.'2 is of size 16x16 and the rest are 
32 x 32 pictures'. Fig. 6 shows the result of reconstructing 
Fig. 1 without averaging. The reconstruction is poor in 
quality and the MNSQR error is 0.0587 which is comparatively 
large. Figs. 7 to 11 show reconstructions with averaging. 
The results are tabulated in Table 2« 
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TABLE 2 


Serial No.' 

Original 

picture 

Reconstructed 

picture 

Size 

MNSQR 

1 

Fig.l 

Fig. 7 

3 2x32 

0.0441 

2 

Fig. 2 

Fig. 8 

16x16 

0.0519 

3 

Fig. 3 

Fig. 9 

32x32 

0.0461 

4 

Fig .4 

Fig. 10 

32x32 

0.0408 

5 

Fig. 5 

Fig.'ll 

32x32 

0.0404 


Table 3 shows the results of increasing the number of 
equations by increasing the number of angles. The original 
picture is Fig. IV 


TABLE 3 


S'.'No'"’ 

No.of 

angles 

Rays per 
angle 

No. of 
equations 

Reconstru- MNSQR 
cted pic. 

1 

13 

32 

416 

Fig. 12 

0';'0322 

2 

15 

32 

480 

Fig. 13 

0.0383 

3 

17 

32 

544 

Fig. '7 

0.0441 

4' 

19 

32 

608 

Fig. 14 

0.0465 

5 

21 

32 

672 

Fig; 15 

0.0510 


The actual C.P.U. time required for one iteration of 
our method is 0,51 seconds which compares favourably with other 


methods’.' 
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The averages considered so far have been 8 point 
averages”*' 4 point averages (Fig. 5.1 and 5.2) were also 
considered. The reconstructions are shown in Figs. 16 and 
17.' The MNSQR errors are 0.0457 and 0.03 94. 

Fig. 18 is the result of doing median filtering at 
each iteration of residual projection'.'. The result is good/ 
especially from the edge preservation point of view. 

* t 

Figs'.' 19 -22 show the results of the various enhance- 
ment schemes";' These enhancement schemes were carried out for 
the reconstruction of Fig. 7. The first three pictures are 
enhancements due to averaging with thresholding/ filtering 
by using the Laplacian filter and removal of negative entries'.' 
The enhanced pictures are poor in quality'.' Only the edge 
preserving smoothing method gives a reasonably good enhanced 
picture as can be seen from Fig. 22. 

In conclusion one can say that if a picture of good 
quality is required/ it is better to increase the number of 
angles than to use an enhancement scheme. This is because 
the successful implementation of any enhancement scheme 
requires an accurate foreknowledge of the mathematical nature 
of the degradation";' The degradation resulting while trying 

to reconstruct a picture with limited projections has not yet 

# 

been mathematically characterized. 
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Fig. 23 -26 show the result of simulations with noise. 
Fig. 23 is the picture of uniformly distributed random noise 
between 0 and 1. The result of reconstructing this is shown 
in Fig. 24.' The reconstructed picture shows smoothing of noise 
due to the process of averrging. This noise is added to the 
original picture and the result is shown in Fig. 25. Fig. 26 
is the reconstruction of the noisy picture. It is to be 
noted, that though it appears as if some lines in the dark 
portion of the picture have been removed/ the fact is that 
the character substituting for these lines (°) has a intensity 
value greater than the characters of the oackground (Z/) # x 
etc'.') . A better overprinting subroutine would surely have 
resulted in a more visually satisfying picture. 
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CHAPTER 6 
CONCLUSION 

A powerful iterative method namely the Residual 
Projection Algorithm has been used to reconstruct pictures 
from their projections. Due consideration has always been 
given to the fact, that for practical reasons, the number 
of projections is always much less than the total npnber 
of pixels. This gave rise to poor reconstructions and the 
problem "was solved by adding physical constraints in the 
form of averaging and median filtering. The resulting 
picture was good in quality 7 

The computational efforts for each of the various 
popular methods have been calculated and it has been shown 
that our method compares favourably in this respect too. 

It requires the same order of computations as ART and is 
free from operations like back projection. Back projection, 

i 

it may be mentioned, is the step that slows down the other- 
wise fast Fourier Reconstruction methods';' An additional 
useful feature of this method is its ability to reconstruct 
in the presence of noise. 

Enhancement schemes have been tried, with very little 
success, to improve the picture quality. The degradation of 
a reconstructed picture has yet to be mathematically 
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characterized. Until this is done, it will be very difficult 
to find a suitable enhancement scheme to improve the quality 
of the reconstructed picture. Till then, the only way of 
obtaining a better picture will be to reconstruct it with an 
increase in the number of angles and hence the number of 
equations".' 

Based on our experiments we suggest any one of the 
two methods, l) Residual projection with averaging, 
preferably the 4 point average of Fig. 5.2, at each itera- 
tion, 2) Residual projection with median filtering at each 
iteration. Both of these methods yield satisfactory 
results'." 

Here, then, we have a method which, because of its 
block iterative nature, is useful in handling large systems, 
is reasonably fast and results in reconstructions which are 
favourably comparable with reconstructions from other 
algorithms In fact, it satisfies all the requirements of 
a ‘good* algorithm. 
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